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Abstract 

The notion of stiffness, which originated in several applications of a different 
nature, has dominated the activities related to the numerical treatment of dif- 
ferential problems for the last fifty years. Contrary to what usually happens in 
Mathematics, its definition has been, for a long time, not formally precise (ac- 
tually, there are too many of them). Again, the needs of applications, especially 
those arising in the construction of robust and general purpose codes, require 
nowadays a formally precise definition. In this paper, we review the evolution 
of such a notion and we also provide a precise definition which encompasses all 
the previous ones. 

Keywords: stiffness, ODE problems, discrete problems, initial value problems, 
boundary value problems, boundary value methods. 



1 Introduction 

Frustra fit per plura quod potest per pau- 
ciora. 

Razor of W. of Ockham, doctor invincibilis. 

The struggle generated by the duality short times-long times is at the heart of 
human culture in almost all its aspects. Here are just a few examples to fix the idea: 



• in historiography: Braudel's distinction among the geographic, social and indi- 
vidual timesjj 
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• in the social sphere: Societies are organized according to three kinds of laws, 
i.e., codes (regulating short term relations), constitutions (regulating medium 
terms relations), and ethical laws (long term rules) often not explicitly stated 
but religiously accepted; 

• in the economy sphere: the laws of this part of human activities are partially 
unknown at the moment. Some models (e.g., the Goodwin model [19]), permits 
us to say, by taking into account only a few variables, that the main evolution is 
periodic in time (and then predictable), although we are experiencing an excess 
of periodicity (chaotic behavior). Nevertheless, some experts claim (see, e.g., 
[T8] ) that the problems in the predictability of the economy are mainly due to 
a sort of gap in passing information from a generation to the next ones, i.e. to 
the conflict between short time and long time behaviors^] 

Considering the importance of this concept, it would have been surprising if the 
duality "short times-long times" did not appear somewhere in Mathematics. As a 
matter of fact, this struggle not only appears in our field but it also has a name: 
stiffness. 

Apart from a few early papers [TQl [IT] , there is a general agreement in placing 
the date of the introduction of such problems in Mathematics to around 1960 |17j . 
They were the necessities of the applications to draw the attention of the mathemat- 
ical community towards such problems, as the name itself testifies: "they have been 
termed stiff since they correspond to tight coupling between the driver and the driven 
components in servo-mechanism" ([12] quoting from |llj). 

Both the number and the type of applications proposing difficult differential prob- 
lems has increased exponentially in the last fifty years. In the early times, the 
problems proposed by applications were essentially initial value problems and, conse- 
quently, the definition of stiffness was clear enough and shared among the few experts, 
as the following three examples evidently show: 

Dl : Systems containing very fast components as well as very slow components 
(Dahlquist [lj). 

D2 : They represent coupled physical systems having components varying with very 
different times scales: that is they are systems having some components varying 
much more rapidly than the others (Liniger [31] . translated from French). 

D3 : A stiff system is one for which \ m ax is enormous so that either the stability 
or the error bound or both can only be assured by unreasonable restrictions 
on h. . .Enormous means enormous relative to the scale which here is i (the 
integration interval). . . (Miranker [32]). 

2 Even Finance makes the distinction between short time and long time traders. 
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The above definitions are rather informal, certainly very far from the precise 
definitions we are accustomed to in Mathematics, but, at least, they agree on a 
crucial point: the relation among stiffness and the appearance of different time-scales 
in the solutions (see also [23]). 

Later on, the necessity to encompass new classes of difficult problems, such as 
Boundary Value Problems, Oscillating Problems, etc., has led either to weaken the 
definition or, more often, to define some consequence of the phenomenon instead 
of defining the phenomenon itself. In Lambert's book [29] five propositions about 
stiffness, each of them capturing some important aspects of it, are given. As matter 
of fact, it has been also stated that no universally accepted definition of stiffness exists 
[36]. 

There are, in the literature, other definitions based on other numerical difficulties, 
such as, for example, large Lipschitz constants or logarithmic norms |37j . or non- 
normality of matrices |23j . Often is not even clear if stiffness refers to particular 
solutions (see, e.g. [23]) or to problems as a whole. 

Sometimes one has the feeling that stiffness is becoming so broad to be nearly 
synonymous of difficult. 

At the moment, even if the old intuitive definition relating stiffness to multiscale 
problems survives in most of the authors, the most successful definition seems to be 
the one based on particular effects of the phenomenon rather than on the phenomenon 
itself, such as, for example, the following almost equivalent items: 



D4 : Stiff equations are equations where certain implicit methods . . . perform better, 
usually tremendous better, than explicit ones [TT] . 

D5 : Stiff equations are problems for which explicit methods don't work |21j . 

D6 : // a numerical method with a finite region of absolute stability, applied to a 
system with any initial condition, is forced to use in a certain interval of in- 
tegration a step length which is excessively small in relation to the smoothness 
of the exact solution in that interval, then the system is said to be stiff in that 
interval [29|. 

As usually happens, describing a phenomenon by means of its effects may not be 
enough to fully characterize the phenomenon itself. For example, saying that fire is 
what produces ash, would oblige firemen to wait for the end of a fire to see if the ash 
has been produced. In the same way, in order to recognize stiffness according to the 
previous definitions, it would be necessary to apply first on^| explicit method and see 
if it works or not. Some authors, probably discouraged by the above defeats in giving 
a rigorous definition, have also affirmed that a rigorous mathematical definition of 
stiffness is not possible [20J. 

It is clear that this situation is unacceptable for at least two reasons: 

3 It is not clear if one is enough: in principle the definition may require to apply all of them. 
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• it is against the tradition of Mathematics, where objects under study have to 
be precisely defined; 

• it is necessary to have the possibility to recognize operatively this class of prob- 
lems, in order to increase the efficiency of the numerical codes to be used in 
applications. 

Concerning the first item, our opinion is that, in order to gain in precision, it would 
be necessary to revise the concept of stability used in Numerical Analysis, which is 
somehow different from the homonym concept used in all the other fields of Mathe- 
matics, where stable are equilibrium points, equilibrium sets, reference solutions, etc., 
but not equations or problems^ (see also [T7j and [30]). 

Concerning the second item, operatively is intended in the sense that the definition 
must be stated in terms of numerically observable quantities such as, for example, 
norms of vectors or matrices. It was believed that, seen from the applicative point 
of view, a formal definition of stiffness would not be strictly necessary: Complete 
formality here is of little value to the scientist or engineer with a real problem to solve 

Nowadays, after the great advance in the quality of numerical co des@ the use- 
fulness of a formal definition is strongly recognised, also from the point of view of 
applications: One of the major difficulties associated with the study of stiff differen- 
tial systems is that a good mathematical definition of the concept of stiffness does not 
exist [6]. 

In this paper, starting from ideas already partially exposed elsewhere [2j HI [26] , we 
will try to unravel the question of the definition of stiffness and show that a precise 
and operative definition of it, which encompasses all the known facets, is possible. 

In order to be as clear as possible, we shall start with the simpler case of initial 
value for a single linear equation and gradually we shall consider more general cases 
and, eventually, we shall synthesize the results. 

2 The asymptotic stability case 

For initial value problems for ODEs, the concept of stability concerns the behavior 
of a generic solution y(t), in the neighborhood of a reference solution y(t), when 
the initial value is perturbed. When the problem is linear and homogeneous, the 
difference, e(t) = y(t) — y(t), satisfies the same equation as y{t). For nonlinear 
problems, one resorts to the linearized problem, described by the variational equation, 
which, essentially, provides valuable information only when y(t) is asymptotically 

4 Only in particular circumstances, for example in the linear case, it is sometimes allowed the 
language abuse: the nonlinear case may contain simultaneously stable and unstable solutions. 

5 A great deal of this improvement is due to the author of the previous sentence. 
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stable. Such a variational equation can be used to generalize to nonlinear problems 
the arguments below which, for sake of simplicity, concerns only the linear case. 

Originally, stiffness was almost always associated with initial value problems 
having asymptotically stable equilibrium points (dissipative problems) (see, e.g., 
Dahlquist p3]). We then start from this case, which is a very special one. Its pecu- 
liarities arise from the following two factsQ 

• it is the most common in applications; 

• there exists a powerful and fundamental theorem, usually called Stability in 
the first approximation Theorem or Poincare-Liapunov Theorem, along with 
its corollary due to Perron^, which allows us to reduce the study of stability of 
critical points, of a very large class of nonlinearities, to the study of the stability 
of the corresponding linearized problems (see, e.g., [9| [271 [35| [38]). 

The former fact explains the pressure of applications for the treatment of such 
problems even before the computer age. The latter one provides, although not always 
explicitly recognized, the mathematical solid bases for the profitable and extensive 
use, in Numerical Analysis, of the linear test equation to study the fixed- /i stability 
of numerical methods. 

We shall consider explicitly the case where the linearized problem is autonomous, 
although the following definitions will take into account the more general case. 

Our starting case will then be that of an initial value problem having an asymp- 
totically stable reference solution, whose representative is, in the scalar case, 

y' = Xy, te[0,T\, ReA < 0, (1) 
2/(0) = V, 

where the reference solution (an equilibrium point, in this case) has been placed at 
the origin. From what is said above, it turns out that it is not by chance that it 
coincides with the famous test equation. 

Remark 2.1 It is worth observing that the above test equation is not less general 
than y' = Xy + g{t), which very often appears in the definitions of stiffness: the only 
difference is the reference solution, which becomes y{t) = J e x ^ s " l g(s)ds, but not 
the topology of solutions around it. This can be easily seen by introducing the new 
variable z(t) = y{t)—y(t) which satisfies exactly equation (QJ) and then, trivially, must 

6 We omit, for simplicity, the other fact which could affect new definitions, i.e., the fact that the 
solutions of the linear equation can be integrated over any large interval because of the equivalence, 
in this case, between asymptotic and exponential stability. 

7 It is interesting to observe that the same theorem is known as the Ostrowsky 's Theorem, in the 
theory of iterative methods. 
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share the same stiffness. Once the solution z(t) of the homogeneous equation has been 
obtained, the solution y(t) is obtained by adding to it y(t) which, in principle, could 
be obtained by means of a quadrature formula. This allows us to conclude that if any 
stiffness is in the problem, this must reside in the homogeneous part of it, i.e., in 
problem (Qp. 



Remark 2.2 We call attention to the interval of integration [0,T], which depends on 
our need for information about the solution, even if the latter exists for all values oft. 
This interval must be considered as datum of the problem. This has been sometimes 
overlooked, thus creating some confusion. 



Having fixed problem ([T]), we now look for a mathematical tool which allows us to 
state formally the intuitive concept, shared by almost all the definitions of stiffness: 
i.e., we look for one or two parameters which tells us if in [0, T] the solution varies 
rapidly or not. This can be done easily by introducing the following two measures for 
the solution of problem ([1]): 



T 



1 11. 







which, in the present case, assume the values: 



1 1 T* 

7c ~ \ReX\T [ ' \Re\\T ~ T ' 

where T* = |i?eA| _1 is the transient time. The two measures k c , 7 c are called condi- 
tioning parameters because they measure the sensitivity of the solution subject to a 
perturbation of the initial conditions in the infinity and in the l\ norm. 

Sometimes, it would be preferable to use a lower value of j c , i.e., 



7c = j^- (3) 

This amounts to consider also the oscillating part of the solution (see also Remark 12.41 
below). 

By looking at Figured! one realizes at once that a rapid variation of the solution 
in [0, T] occurs when k c 3> 7 C . It follows then that the parameter 



which is the ratio between the two characteristic times of the problem, is more sig- 
nificant. Consequently, the definition of stiffness follows now trivially: 



Definition 2.1 The initial value problem (QP is stiff if cr c 3> 1. 
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Figure 1: Solutions and values of k c and 7 C in the cases A = —.2 (left plot) and A = —2 
(right plot). 

The parameter a c is called stiffness ratio. 

Remark 2.3 The width of the integration interval T plays a fundamental role in the 
definition. This is an important point: some authors, in fact, believe that stiffness 
should concern equations; some others believe that stiffness should concern problems, 
i.e., equations and data. We believe that both statements are partially correct: stiffness 
concerns equations, integration time, and a set of initial data (not a specific one of 
them). Since this point is more important in the non scalar case, it will be discussed 
in more detail later. 

Remark 2.4 When 7 C is defined according to the definition of stiffness contin- 
ues to be also meaningful in the case ReX = 0, i.e., when the critical point is only 
marginally stable. In fact, let 



and the definition encompasses also the case of oscillating stiffness introduced by some 
authors (e.g., JEM/). Once again the stiffness is the ratio of two times. If information 
about the solution on the smaller time scale is needed, an adequately small stepsize 
should be used. It is worth noting that high oscillating systems (with respect to T) 
fall in the class of problems for which explicit methods do not work, and then are stiff 
according to definitions D4-D6. 




Then 



T 



°c = 2tv—, 
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When X = 0, then k c = 7 C = o~ c = 1. 

In the case ReX > (i.e., the case of an unstable critical point) , both parameters k c 
and 7 C grow exponentially with time. This implies that small variations in the initial 
conditions will imply exponentially large variations in the solutions, both pointwise 
and on average: i.e., the problem is ill conditioned. 

Of course, the case ReX = considered above cannot be considered as representa- 
tive of more difficult nonlinear equations, since linearization is in general not allowed 
in such a case. 

The linearization is not the only way to study nonlinear differential (or difference) 
equations. The so called Liapunov second method can be used as well (see, e.g., [22l [27] 
138]). It has been used, in connection with stiffness in (5J [131 El [El EEl Ej, although 
™t : always explicitly named! Anyway, no matter how the asymptotic stability of a 
reference solution is detected, the parameters ([2]) and Definition 12.11 continue to be 
valid. Later on, the problem of effectively estimating such parameters will also be 
discussed. 



2.1 The discrete case 

Before passing to the non scalar case, let us now consider the discrete case, where 
some interesting additional considerations can be made. Here, almost all we have 
said for the continuous case can be repeated. The first approximation theorem can 
be stated almost in the same terms as in the continuous case (see e.g. |28j). 

Let the interval [0, T] be partitioned into N subintervals of length h n > 0, thus 
defining the mesh points: t n = YTj=i hj, n = 0,1, N. 

The linearized autonomous problem is now: 



y n+ i = n n y n , n = 0, . . . , JV - 1, Vo = V, (5) 

where the {/^ n } are complex parameters. The conditioning parameters for (jSJ), along 
with the stiffness ratio, are defined as: 



«d = rr max \vi\i id = A^y^ max (N' l^-il)> o- d = —. (6) 

\V\ i=o,-,N \V\ T ~[ Id 

This permits us to define the notion of well representation of a continuous problem 
by means of a discrete one. 



8 Often, it appears under the name of one-sided Lipschitz condition. 
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method 




condition 


Explicit Euler 


1 + hX 




l + h\\<l 


Implicit Euler 


1 

l-h\ 




i 

1 1-hX \ 


< 1 


Trapezoidal Rule 


l+hX/2 
l-hX/2 




l+hX/2 
l-hX/2 


< 1 



Table 1: Condition (JTj) for some popular methods. 



Definition 2.2 The problem (T7]j is well represented by (TJJ) if 

k c « fc d , (7) 
7c ~ 7d- (8) 

In the case of a constant mesh-size h, \x n = \x and it easily follows that the 
condition (J7|) requires < 1. It is not difficult to recognize the usual A-stability 
conditions for one-step methods (see Tabled]). Furthermore, it is easily recognized 
that the request that condition (j7|) holds uniformly with respect to hX G C~ implies 
that the numerical method producing (jSj) must be implicit. 

What does condition flS} require more? Of course, it measures how faithfully the 
integral J Q T \ y(t)\ dt is approximated by the quadrature formula Xlili ^» max( | ?/j | , | | ) , 
thus giving a sort of global information about the behavior of the method producing 
the approximations {t/i}- One of the most efficient global strategies for changing the 
stepsize is based on monitoring this parameter [3j HJ [7J [8] [33j [M] . In addition to this, 
when finite precision arithmetic is used, then an interesting property of the parameter 
7d occurs [26]: if it is smaller than a suitably small threshold, this suggests that we 
are doing useless computations, since the machine precision has already been reached. 

2.2 The non scalar case 

In this case, the linearized problem to be considered is 

y' = Ay, te[0,T\, y(0)=v, (9) 

with A G ]R mxm and having all its eigenvalues with negative real part. It is clear 
from what was said in the scalar case that, denoting by $(£) = e At the fundamental 
matrix of the above equation, the straightforward generalization of the definition of 
the conditioning parameters fl2]) would lead to: 

k c = max ||$(t)||, 7c = 1 / \Mt)\\dt, a c = ^. (10) 

Indeed, these straight definitions work most of the time, as is confirmed by the fol- 
lowing example, although, as we shall explain soon, not always. 
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Figure 2: Estimated stiffness ratio of Van der Pol's problem (TTTj) . 

Example 2.1 Let us consider the well-known Van der Pol's problem, 

y'i = 2/2, 

y' 2 = -y 1 + f iy 2 {l-y 2 1 ), *e[0,2//], (11) 
2/(0) = (2, Of, 

whose solution approaches a limit cycle of period T rs 2/i. i£ a/so very well-known 
that, the larger the parameter \i, the more difficult the problem is. In Figure [H we 
plot the parameter a c {p) (as defined in [Tfy) ) for ft ranging from to 10 3 . Clearly, 
stiffness increases with fi. 

Even though (TTDT) works for this problem, this is not true in general. The problem is 
that the definition of stiffness as the ratio of two quantities may require a lower bound 
for the denominator. While the definition of k c remains unchanged, the definition of 
7 C is more entangled. Actually, we need two different estimates of such a parameter: 

• an upper bound, to be used for estimating the conditioning of the problem in 
li norm; 

• a lower bound, to be used in defining a c and, then, the stiffness. 

In the definition given in [21 Sj, this distinction was not made, even though the 
definition was (qualitatively) completed by adding 

"/or at least one of the modes". (12) 
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We shall be more precise in a moment. In the meanwhile, it is interesting to note 
that the clarification contained in (fT2l is already in one of the two definitions given 
by Miranker |32j : 

A system of differential equations is said to be stiff on the interval (0, i) if there 
exists a solution of that system a component of which has a variation on that interval 
which is large compared to -=, 

where it should be stressed that the definition considers equations and not problems: 
this implies that the existence of largely variable components may appear for at least 
one choice of the initial conditions, not necessary for a specific one. 

Later on, the definition was modified so as to translate into formulas the above 
quoted sentence ( TT2l . The following definitions were then given (see, e.g., [26]): 



k c (T,T]) = -rr-rr max \\y(t)\\, k c {T) = max/? c (T, rj), 



\V\ 



0<t<T 



T 



(13) 



lc(T,v) = Tfrr^ / \\y(t)\\dt } 7 C (T) = max7 c (T,7?). 
T \\v\\ Jo v 



and 



/rri\ K C {T, Tj) 

a c (T) = max 14 

The only major change regards the definition of a c . Let us be more clear on this 
point with an example, since it leads to a controversial question in the literature: i.e., 
the dependence of stiffness from the initial condition. Let A = diag(Ai, A 2 , . . . , A m ) 
with Aj < and |Ai| > |A 2 | > . . . > |A m |. The solution of problem ([9]) is y(t) = e At rj. 

If a c is defined according to ({TO]) , it turns out that ||e A '|| = e Xmt and, then, 
7 C (T) « T\x m \ ■ ^' nowever 5 we take i] = (1, 0, . . . , 0) T , then y(t) = e Xlt and J C (T) 
becomes 7 C (T) w j^j- Of course, by changing the initial point, one may activate 

each one of the modes, i.e. the functions e Xit on the diagonal of the matrix e At , leaving 
silent the others. This is the reason for specifying, in the older definition, the quoted 
sentence (1121) . The new definition (Tl4|) . which essentially poses as the denominator of 
the ratio a c the smallest value among the possible values of 7 C (T, rj), is more compact 
and complies with the needs of people working on the construction of codes, who 
like more operative definitions. For the previous diagonal example, we have that k c 
continues to be equal to 1, while 7 C (T) = ^p^j- 

Having got the new definition (Ti"4"l) of cr c (T), the definition of stiffness continues 
to be given by Definition 12.11 given in the scalar case, i.e., the problem (Q is stiff if 

<7 C (T)»1. 

How does this definition reconcile with the most used definition of stiffness for 
the linear case, which considers the "smallest" eigenvalue A m as well? The answer is 
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already in Miranker's definition D3. In fact, usually the integration interval is chosen 
large enough to provide complete information on the behavior of the solution. In this 
case, until the slowest mode has decayed enough, i.e. T = l/|A m |, which implies 



Or 



T 



Ai 



A, 



(15) 



which, when much larger than 1, coincides with the most common definition of stiff- 
ness in the linear case. However, let us insist on saying that if the interval of integra- 



tion is much smaller than l/|A m |, the 'problem may be not stiff even if 



> 1. 



The controversy about the dependence of the definition of stiffness on the initial 
data is better understood by considering the following equation given in [29], PP- 217- 
2181: 



dt\y 2 J V -1-999 °- 999 ) \V2 ) \ 0.999(sint-cost) 
whose general solution is 



(S)-^(i) + ^(ii.) + («0- 

The initial condition y(0) = (2, 3) T requires c 2 = and, then, the slowest mode is 
not activated: the solution rapidly reaches the reference solution. If this information 
was known beforehand, one could, in principle, choose the interval of integration T 
much smaller than q-^-. This, however, does not take into account the fact that 
the computer uses finite precision arithmetic, which may not represent exactly the 
initial condition rj. To be more precise, let us point out that the slowest mode is 
not activated only if the initial condition is on the line 1/2(0) — 2/i(0) — 1 = 0. Any 
irrational value of yi(0) will not be well represented on the computer. This is enough 
to activate the silent mode. Of course, if one is sure that the long term contribution 
to the solution obtained on the computer is due to this kind of error, a small value of 
T can always be used. But it is rare that this information is known in advance. For 
this reason, we consider the problem to be stiff, since we believe that the definition of 
stiffness cannot distinguish, for example, between rational and irrational values of the 
initial conditions. Put differently, initial conditions are like a fuse that may activate 
stiffness. 

We conclude this section by providing a few examples, which show that Defini- 
tion 12.11 when a c is defined according to (|T4l) . is able to adequately describe the 
stiffness of nonlinear and/or non autonomous problems as well. 
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Figure 3: Estimated stiffness ratio of Robertson's problem (f!6l) . 

Example 2.2 Let us consider the well-known Robertson's problem: 
y[ = -M yi + I0 4 y 2 y 3 , 

y' 2 = M yi -10 4 y 2 y 3 -3-10 7 yl t e [0,T], (16) 

y' 3 = 3-io 7 y 2 2 , 

y(0) = (1, 0, 0) T . 

Its stiffness ratio with respect to the length T of the integration interval, obtained 
through the linearized problem and considering a perturbation of the initial condition 
of the form (0, e, —s) T , is plotted in Figured As it is well-known, the figure confirms 
that for this problem stiffness increases with T . 

Example 2.3 Let us consider the so-called Kreiss problem l21\ p. 542], a linear and 
non autonomous problem: 

y' = A(t)y, te[0,47r], y(0) fixed, (17) 

where 

A{t) = Q T (t)k £ Q(t), (18) 

and 

/ cost sint \ . ( —\ \ 

«W=(-smi cos* )• A '=( - £ -J' (9) 

Its stiffness ratio with respect to the small positive parameter e, obtained by consider- 
ing a perturbation of the initial condition of the form (— e, 1) T , is plotted in Figure^ 
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Example 2.4 Let us consider the following linear and non autonomous problem, a 
modification of problem [Fty, that we call "modified Kreiss problem" $ 



where 



and 



y' = A(t)y, te[0,47r], y(0) fixed, (20) 
A(t) = Q-\t)p- l KPQ £ (t), (21) 



P=( } j), QM=(jnt Jnt), A«=f 1 _ £ -lj. (22) 

Its stiffness ratio with respect to the small positive parameter e, obtained by consider- 
ing a perturbation of the initial condition of the form (— e, 1) T , is shown in Figured 
Also in this case the stiffness of the the problem behaves as e~ x , as e tends to 0. 



Remark 2.5 It is worth mentioning that, in the examples considered above, we nu- 
merically found that 



max 



' This problem has been suggested by J.I. Montijano. 
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Figure 5: Estimated stiffness ratio of the modified Kreiss problem ( f20l) - ([22l) . 

is obtained by considering an initial condition rj in the direction of the eigenvector 
of the Jacobian matrix (computed for t « t ) associated to the dominant eigenvalue. 
We note that, for an autonomous linear problem, if A is diagonalizable, this choice 
activates the mode associated with Ax, i.e., the eigenvalue of maximum modulus of A. 

2.3 The non scalar discrete case 

As for the scalar case, what we said for the continuous problems can be repeated, 
mutatis mutandis, for the discrete ones. For brevity, we shall skip here the details for 
this case, also because they can be deduced from those described in the more general 
case discussed in the next section. 

3 Boundary Value Problems (BVPs) 

The literature about BVPs is far less abundant than that about IVPs, both in the 
continuous and in the discrete case. While there are countless books on the latter 
subject presenting it from many points of view (e.g., stability of motion, dynamical 
systems, bifurcation theory, etc.), there are many less books about the former. More 
importantly, the subject is usually presented as a by product of the theory of IVPs. 
This is not necessarily the best way to look at the question, even though many 
important results can be obtained this way. However, it may sometimes be more 
useful to look at the subject the other way around. Actually, the question is that 
IVPs are naturally a subclass of BVPs. Let us informally clarify this point without 
many technical details which can be found, for example, in [I]. 
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IVPs transmit the initial information "from left to right" . Well conditioned IVPs 
are those for which the initial value, along with the possible initial errors, decay 
moving from left to right. FVPs (Final Value problems) are those transmitting in- 
formation "from right to left" and, of course, well conditioning should hold when the 
time, or the corresponding independent variable, varies towards — oo. More precisely, 
considering the scalar test equation ([!]), the asymptotically stability for IVPs and 
FVPs requires ReX < and ReX > 0, respectively. BVPs transmit information both 
ways. Consequently, they cannot be scalar problems but vectorial of dimension at 
least two. We need then to refer to the test equation (J9]). It can be affirmed that a 
well conditioned linear BVP needs to have eigenvalues with both negative and posi- 
tive real parts (dichotomy, see, e.g., [UH]). More precisely: the number of eigenvalues 
with negative real part has to match the amount of information transmitted "from 
left to right", and the number of eigenvalues with positive real part has to match 
the amount of information traveling "from right to left". For brevity, we shall call 
the above statement continuous matching rule. Of course, if there are no final con- 
ditions, then the problem becomes an IVP and, as we have seen, in order to be well 
conditioned, it must have all the eigenvalues with negative real part. In other words, 
the generalization of the case of asymptotically stable IVPs is the class of well condi- 
tioned BVPs because both satisfy the continuous matching rule. This is exactly what 
we shall assume hereafter. 

Similar considerations apply to the discrete problems, where the role of the imagi- 
nary axis is played by the unit circumference in the complex plane. It is not surprising 
that a numerical method will well represent a continuous autonomous linear BVP if 
the corresponding matrix has as many eigenvalues inside the unit circle as the number 
of initial conditions and as many eigenvalues outside the unit circle as the number of 
final conditions (discrete matching rule). 

Remark 3.1 The idea that IVPs are a subset of BVPs is at the root of the class 
of methods called Boundary Value Methods (BVMs) which permits us, thanks to 
the discrete matching rule, to define high order and perfectly A-stable methods (i.e., 
methods having the imaginary axis separating the stable and unstable domains), which 
overcome the Dahlquist's barriers, and are able to solve both IVPs and BVPs (see, 

e.g., W- 

Remark 3.2 From this point of view, the popular shooting method, consisting of 
transforming a BVP into an IVP and then applying a good method designed for 
IVPs, does not appear to be such a good idea. As matter of fact, even a very well 
conditioned linear BVP, i.e. one which satisfies the continuous matching rule, will be 
transformed in a badly conditioned IVP, since the matrix of the continuous IVP shall, 
of course, contain eigenvalues with positive real part. This will prevent the discrete 
matching rule to hold. 
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3.1 Stiffness for BVPs 



Coming back to our main question, stiffness for BVPs is now defined by generalizing 
the idea already discussed in the previous sections. 

As in the previous cases, we shall refer to linear problems, but the definitions 
will also be applicable to nonlinear problems as well. Moreover, according to what is 
stated above, we shall only consider the case where the problems are well conditioned 
(for the case of ill conditioned problems, the arguments are slightly more entangled, 
see e.g. [7]). Then, let us consider the linear and non autonomous BVP: 

y' = A(t)y, t e [0, T], B y(0) + B x y{T) = V , (23) 
where y(t), i)6l m and A(t),B Q , B x e R mxm . The solution of the problem (EHj) is 

y(t) = $(t)Q-\ 



where $(t) is the fundamental matrix of the problem such that $(0) = /, and Q = 
B a + Bb<&(T), which has to be nonsingular, in order for (!23j) to be solvable!^! 

As in the continuous IVP case, the conditioning parameters are defined (see ([131) ) 

as: 



K c (T,ri) = |^ max ||y(t)||, k c (T) = max k c (T, 77), 

(24) 

1 f T 

lc{T,rj) = —-r / \\y{t)\\dt, 7 C (T) = max7 c (T,r]). 
T \\V\\ Jo 

Consequently, the stiffness ratio is defined as (see (fT4|) ): 

a c {T) = max , 
»» lc[T,V) 

and the problem is stiff if cr c (T) ^> 1. Moreover, upper bounds of k c {T) and 7 C (T) 
are respectively given by: 

k c (T) < max \m)Q-% 7c (T) < ^ F W^Q-^dt. (25) 

0<t<T 1 J Q 

Thus, the previous definitions naturally extend to BVPs the results stated for 
IVPs. In a similar way, when considering the discrete approximation of f)23p . for 
the sake of brevity provided by a suitable one-step method over a partition it of the 



' Observe that, in the case of IVPs, Bq = I and B\ = O, so that Q = I. 
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interval [0, T], with subintervals of length hi, i = 1, . . . , N, the discrete problem will 
be given by 

y n+1 = R n y n , n = 0, . . . , JV - 1, B y + B x y N = 77, (26) 
whose solution is given by 

(n-1 \ N-1 
J] ^ QnV, Qn = B + Bt J] 
i=0 / i=0 

The corresponding discrete conditioning parameters are then defined by: 
K d (ir, r]) = T|— if max lbn||, Kd(7r) = maxK d (7r,^), 

m 0<n<N r] 



(27) 



7d(v,r/) = — — y"/iimax(||yj||, Tdl 71 ") = max7 d (vr, 77), 

1 ?7 Z ' 7) 

II II 



and 

cr d (vr) = max — -. 

According to Definition [2^2], we say that the discrete proble ( T26|) ioe// represents 
the continuous problem fl23|) if 

KdW««c(T), 7 dW «7c(T). (28) 

Remark 3.3 It is worth mentioning that innovative mesh- selection strategies for the 
efficient numerical solution of stiff BVPs have been defined by requiring the match 
m> (^e, e.g., |3 » ' 



3.2 Singular Perturbation Problems 

The numerical solution of singular perturbation problems can be very difficult because 
they can have solutions with very narrow regions of rapid variation characterized by 
boundary layers, shocks, and interior layers. Usually, the equations depend on a small 
parameter, say e, and the problems become more difficult as e tends to 0. It is not 
always clear, however, how the width of the region of rapid variation is related to 
the parameter e. By computing the stiffness ratio o~ c (T), we observe that singularly 
perturbed problems are stiff problems. Moreover, as the following examples show, 
the parameter c c (T) provides us also with information about the width of the region 
of rapid variation. 



It is both defined by the used method and by the considered mesh. 
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Figure 6: Problem (E9l). e= 1(T 8 . 

The examples are formulated as second order equations: of course, they have to 
be transformed into corresponding first order systems, in order to apply the results 
of the previous statements. 

Example 3.1 Let us consider the linear singularly perturbed problem: 

£y " + ty' = -en 2 cos(?rt) - nt sin(Trf), y(-l) = -2, y(l) = 0, (29) 

whose solution has, for < e <C 1, a turning point at t = (see Figure^ . The exact 
solution is y{t) = cos(nt) + exp((t — \)/ yfe) + exp(— (t + 1)/ yfe). 

In Figure [?| we plot an estimate of the stiffness ratio obtained by considering two 
different perturbations of the boundary conditions of the form (1, 0) T and (0, 1) T . The 
parameter e varies from 10 _1 to 10~ 14 . We see that the ( estimated) stiffness parameter 
grows like Ve^ 1 - 



Example 3.2 Let us consider the following nonlinear problem: 

7T / Tit \ 

ey" + exp(y)y' - - sin I- J exp(2y) = 0, y(0) = 0, y(l) = 0. (30) 

This problem has a boundary layer at t = (see Figure^). In Figured we plot an 
estimate of the stiffness ratio obtained by considering two different perturbations of 
the boundary conditions of the form (1, 0) T and (0, 1) T . The parameter e varies from 
1 to 10~ 8 . We see that the (estimated) stiffness parameter grows like e~ x , as e tends 
to 0. 
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Figure 7: Estimated stiffness ratio of problem ( 1291) . 

Example 3.3 Let us consider the nonlinear Troesch problem: 

y" = tiBmh(jjLy), y(Q) = 0, y(l) = 1. (31) 

This problem has a boundary layer near t = 1 (see Figure [70j) . In Figure [771 we p/ot 
i/ie estimate of the stiffness ratio obtained by considering two different perturbations of 
the boundary conditions of the form (1, 0) T and (0, 1) T . The parameter fi is increased 
from 1 to 50 and, as expected, the stiffness ratio increases as well: for \x = 50, it 
reaches the value 1.74 ■ 10 12 . 
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